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Abstract 

The  goal  of  this  study  is  to  assess  models  for  Deep  Convection  with  special  emphasis  on  their  use  in  coarse 
resolution  ocean  general  circulation  models.  A  model  for  deep  convection  must  contain  both  vertical 
transport  and  lateral  advection  by  mesoscale  eddies  generated  by  baroclinic  instabilities.  The  first  process 
operates  mostly  in  the  initial  phases  while  the  second  dominates  the  final  stages.  Here,  the  emphasis  is  on 
models  for  vertical  mixing.  When  mesoscales  are  not  resolved,  they  are  treated  with  the  Gent  and 
McWilliams  parameterization.  The  model  results  are  tested  against  the  measurements  of  Lavender,  Davis 
and  Owens,  2002  (LDO)  in  the  Labrador  Sea.  Specifically,  we  shall  inquire  whether  the  models  are  able  to 
reproduce  the  region  of  “ deepest  convection  ”  which  we  shall  refer  to  as  DC  (mixed  layer  depths  800-1300 
m).  The  region  where  it  was  measured  by  Lavender  et  al.  (2002)  will  be  referred  to  as  the  LDO  region.  The 
main  results  of  this  study  can  be  summarized  as  follows. 

(1)  3°  x  3°  resolution .  A  GFDL-type  OGCM  with  the  GISS  vertical  mixing  model  predicts  DC  in  the  LDO 
region  where  the  vertical  heat  diffusivity  is  found  to  be  10  m2  s“l,  a  value  that  is  quite  close  to  the  one 
suggested  by  heuristic  studies.  No  parameter  was  changed  from  the  original  GISS  model.  However,  the 
GISS  model  also  predicts  some  DC  in  a  region  to  the  east  of  the  LDO  region. 

(2)  3°  x  3°  resolution .  A  GFDL-type  OGCM  with  the  KPP  model  (everything  else  being  the  same)  does  not 
predict  DC  in  the  LDO  region  where  the  vertical  heat  diffusivity  is  found  to  be  0.5  x  10“4  m2  s~!  which  is 
the  background  value.  The  KPP  model  yields  DC  only  to  the  east  of  the  LDO  region. 

(3)  1°  x  1°  resolution.  In  this  case,  a  MY2.5  mixing  scheme  predicts  DC  in  the  LDO  region.  However,  it  also 
predicts  DC  to  the  west,  north  and  south  of  it,  where  it  is  not  observed.  The  behavior  of  the  KPP  and 
MY  models  are  somewhat  anti-symmetric.  The  MY  models  yield  too  low  a  mixing  in  stably  stratified 
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flows  since  they  predict  a  critical  Richardson  number  Ri( cr)  =  0.19  which  is  five  times  smaller  than  the 
value  Ri( cr)  =  0(1)  needed  to  obtain  realistic  ML  depths.  However,  as  discussed  above,  in  unstable 
stratifications  the  MY  models  yield  better  results.  On  the  other  hand,  the  KPP  model,  which  was  mo¬ 
tivated  primarily  by  the  need  to  overcome  the  MY  “too  low  mixing”  in  stable  stratification,  yields  at 
coarse  resolution,  no  DC  in  the  LDO  region.  In  this  respect,  the  GISS  model,  yields  both  a  correct 
Ri(ct)  =  0(1)  in  stable  stratification  and  correct  results  in  the  unstable  configuration  in  the  LDO 
region. 

(4)  1/3°  x  1/3°  resolution.  In  this  case,  KPP  predicts  mixed  layer  depths  up  to  1.7  km  inside  the  LDO  region 
where  at  coarse  resolution  none  existed.  However,  the  model  still  produces  DC  at  locations  outside  the 
LDO  region  where  it  is  not  observed.  However,  since  these  regions  are  intermingled  with  very  shallow 
mixed  layer  depths,  the  resulting  mean  mixed  layer  depths  turn  out  to  be  less  than  800  m  almost  every¬ 
where  outside  the  LDO  region. 

(5)  1/12°  x  1/12°  resolution.  In  this  case,  KPP  predicts  mixed  layer  depths  up  to  3  km  both  inside  and  out¬ 
side  the  LDO  region.  These  regions  are,  here  too,  intermingled  with  very  shallow  mixed  layer  depths 
with  resulting  mean  mixed  depths  greater  than  800  m  both  inside  and  outside  the  LDO  region. 

In  conclusion,  as  for  a  model  for  deep  convection  to  be  used  in  coarse  resolution,  these  results  indicate 
that  the  GISS  mixing  model  fares  well  with  observations  in  both  stable  and  unstable  stratifications  but 
overestimates  its  geographical  extent.  This  leads  to  the  problem  of  future  improvements  of  the  model.  It 
must  be  generalized  to  include  the  following  physically  important  features:  (a)  rotation  that  becomes  im¬ 
portant  in  the  later  phases  of  deep  convection  when  it  acts  to  slow  down  the  rate  of  mixed  layer  deepening, 
(b)  non-locality,  in  particular  skewness  which  is  large  (negative)  in  the  initial  phases  of  deep  convection  and 
becomes  small  in  the  final  stages,  and  finally,  (c)  a  new  model  to  treat  lateral  advection  by  baroclinic  eddies 
that  in  the  final  stages  of  deep  convection  dominates  over  vertical  transport. 

©  2003  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Earth’s  atmosphere  interacts  with  the  ocean  in  two  ways.  Wind  stresses  cause  strong  mixing  in 
the  first  100  m  of  the  ocean  (mixed  layer,  ML)  but  hardly  affect  water  masses  below  the  ML 
where  lies  the  largest  portion  of  the  ocean  characterized  by  stable  stratification  and  thus  weak 
mixing.  To  the  effect  of  communicating  with  the  deep  ocean,  e.g.,  in  the  process  of  absorbing 
atmospheric  C02,  stable  stratification  acts  like  a  rigid  lid  that  insulates  the  strongly  mixed  ML 
from  the  weakly  mixed  deep  ocean.  If  such  a  configuration  always  prevailed,  the  deep  ocean 
would  be  shielded  from  climatic  events;  deep  waters  would  be  dynamically  decoupled  from  surface 
phenomena  and  the  ocean  deep  currents  would  be  considerably  weaker  than  what  is  observed. 
Differential  solar  heating  between  low  and  high  latitudes  would  not  result  in  a  poleward  flow  of 
warm  waters  and  the  Atlantic  would  look  more  like  the  Pacific  where  there  are  no  deep  convective 
regions  (Weaver  et  al.,  1999).  Earth’s  climate  would  be  quite  different  from  what  is  observed 
today. 

Deep  Convection  is  the  only  process  through  which  surface  phenomena  such  as  buoyancy 
losses,  brine  rejection,  etc.,  pierce  the  lid  of  strong  stable  stratification  that  characterizes  the 
thermocline.  Ultimately,  this  leads  to  the  formation  of  deep  waters  (Chu  and  Gascard,  1991; 
Maxworthy,  1997;  Marshall  and  Schott,  1999;  The  Labrador  Sea  Deep  Convection  Experiment, 
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2002).  Open  ocean  deep  convection  occurs  in  a  small  number  of  locations,  Labrador,  Greenland, 
Weddell  and  Western  Mediterranean  Seas,  and  yet  it  is  one  of  the  ocean’s  major  features  since  it 
represents  the  initial  stage  of  the  global-scale  ventilation  loops  of  the  world  ocean  (Sander  et  al., 
1995).  In  fact,  it  is  the  dominant  mechanism  responsible  for  the  production  of  North  Atlantic 
Deep  Water  (NADW)  and  of  the  Antarctic  Bottom  Water  (AABW,  Weaver  et  al.,  1999).  Both 
NADW  and  AABW  play  a  major  role  in  earth’s  climate. 

Loss  of  surface  buoyancy  and/or  brine  rejection  lead  to  a  top-heavy,  unstable  configuration 
which  acts  as  precursor  of  turbulent  motion  that  ultimately  leads  to  deep  convection.  The  latter 
upwells  warmer  waters  that  can  melt  ice  and  reduce  the  albedo  resulting  in  a  “negative  feedback” 
that  affects  climate  (Killworth,  1983).  Regrettably,  however,  deep  conVection  is  still  poorly 
modeled  in  coarse  resolution  OGCM  (Killworth,  1989).  While  laboratory  and  numerical  simu¬ 
lations  (Sander  et  al.,  1995;  Denbo  and  Skyllingstad,  1996;  Skyllingstad  et  al.,  1996;  Maxworthy, 
1997;  Marshall  and  Schott,  1999)  have  brought  to  light  several  key  features  of  convective  pro¬ 
cesses,  the  translation  of  this  information  into  a  reliable  model  for  coarse  resolution  OGCMs  has 
not  yet  been  achieved  but,  given  the  complexity  of  the  phenomenon,  this  is  hardly  surprising. 
Before  we  test  models  for  deep  convection,  it  is  important  to  discuss  some  of  its  key  features: 

(1)  Deep  convection  is  a  highly  turbulent  process.  This  is  exhibited  by  the  large  vertical  diffu- 
sivities  Kv: 

Kv  ~  /w  ~  (1-10)105  cm2s~'  (la) 

where  we  have  used  /  ~  (1-2)  km  and  w  ~  (1-5)  cms'1,  as  discussed  by  Marshall  and  Schott 
(1999).  Equivalently,  one  may  consider  the  large  value  of  the  Reynolds  number: 

Re~Kv/v~  107  (lb) 

where  v  ~  10-2  cm2  s_I  is  the  kinematic  viscosity  of  seawater.  Thus,  a  high-Re  turbulence  model 
is  needed  to  describe  deep  convective  processes. 

(2)  Deep  convection  is  affected  by  rotation.  Consider  the  characteristic  length  scale  (Golystin, 
1980), 

/(rot)  =  (5,//3)1/2  -  (0.15-0.56)  km  (lc) 

where  B »  is  the  surface  buoyancy  and  /  is  the  Coriolis  parameter.  The  numerical  values  in  (lc) 
correspond  to  the  Greenland  and  Mediterranean  Seas  (Marshall  and  Schott,  1999,  Table  3.4.1). 
Contrary  to  the  atmospheric  case  whose  -/(rot)  is  much  larger  than  the  height  «1  km  of  the 
planetary  boundary  layer,  in  the  ocean  case  the  reverse  is  often  true,  namely  /(rot)  may  be 
considerably  smaller  than  the  ocean  depth  H  that,  for  the  two  cases  just  cited,  is  1.5  and  1.8  km. 
This  yields  small  Rossby  numbers  Ro  =  /(rot )/H  =  0. 1-0.3,  an  indication  of  the  importance  of 
rotation. 

A  turbulence  model  must  be  able  to  incorporate  rotational  effects  and  more  specifically,  it  must 
be  able  to  reproduce  key  features  like  the  Golystin’s  scale  (lc).  Rotation  enters  the  turbulence 
equations  not  only  through  the  linear  Coriolis  term  in  the  dynamic  equations  for  the  turbulent 
velocity  but,  more  important,  it  affects  the  very  structure  of  the  non-linear  interactions  that  are  at 
the  heart  of  turbulence.  In  the  presence  of  rotation,  velocity  components  with  different  vectors  are 
rotated  by  the  Coriolis  force  around  different  axes  that  coincide  with  the  directions  of  the  cor¬ 
responding  wave-vectors.  Thus,  the  energy  cascade  from  large  to  small  eddies  is  inhibited.  An 
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inertial  range  is  still  present  but  only  for  wave  numbers  larger  than  A:(rot)  where  the  latter  is  the 
inverse  of  Eq.  (lc).  For  wave  numbers  k  <  k{ rot),  the  spectrum  is  no  longer  of  the  Kolmogorov 
type.  That  is,  one  has  two  regimes: 

£  >  ^(rot)  :  E(k)  ~  (eQ)l,2k~2,  £  <  f(rot)  :  E{k)  ~  e2/3k~5/3  (Id) 

where  £  is  the  rate  of  energy  dissipation.  Integrating  the  two  spectra,  one  derives  that  the  cor¬ 
responding  velocities  with  and  without  rotation  are  derived  to  be 

w(rot)  ~  [£/£(rot)}l/2(Bff)l/2,  w„  ~  (le) 

where  we  have  taken  the  dissipation  rate  e  equal  to  the  surface  buoyancy  flux  Bt.  Values  of  w(rot) 
and  w0  for  the  Mediterranean,  Labrador  and  Greenland  Seas  can  be  found  in  Table  3.4.1  of 
Marshall  and  Schott  (1999).  A  turbulence  model  capable  of  predicting  the  two  regimes  of  the 
energy  spectrum  (Id)  has  recently  been  constructed  and  its  implications  tested  against  direct 
numerical  simulations  (Canuto  and  Dubovikov,  1997). 

(3)  Deep  convection  cannot  be  fully  represented  by  a  purely  local  model.  If  a,  1  -  er  are  the  areas 
occupied  by  the  up/downdrafts  (plumes),  one  has  (Moeng  and  Rotunno,  1990): 

o-  =  1  /2(1  -  D),  D  =  Sw(l  +  Sl)~l/2,  Sw  =  mP/w?'2  (2a) 

where  Sw  is  the  skewness  of  the  velocity  field.  Labrador  Sea  data  (Lavender  et  al.,  2002;  Herbaut 
and  Marshall,  2002)  show  that  in  the  initial  stages  of  convective  deepening,  Sw  is  large  and 
negative  while  in  the  final  stages  it  is  small:  skewness  subsides  with  time  (for  numerical  values  of 
Sw,  see  Table  1  of  Lavender  et  al.,  2002).  This  means  that  updrafts  subside  with  time  and  that  the 
downward  vertical  velocities  occur  more  frequently  than  upward  velocities.  By  using  a  local  model 
with  Sw  =  0,  or  equivalently: 

,  =  1/2  (2b) 

one  assumes  that  downdrafts  and  updrafts  occupy  the  same  fractional  areas,  that  is,  a  local  model 
represents  only  the  final  stages  of  deep  convection.  On  the  other  hand,  the  OPPS  model  (Ocean 
Penetrative  Plume  Scheme)  of  Paluszkiewicz  et  al.  (1994),  Paluszkiewicz  and  Romea  (1997)  rep¬ 
resents  the  initial  stages  of  convective  processes.  In  fact,  consider  the  two  assumptions  underlying 
the  OPPS;  (a)  an  isolated  plume.  This  means  that  the  fractional  area  occupied  by  the  plume  must 
be  small, 

1  —  a  <C  1  (2°) 

From  (2a)  it  follows  that  Sw  must  be  large  and  negative  and  thus  a  plume  model  applies  only  to  the 
initial  stages  of  deep  convection;  (b)  a  quiescent  environment  in  which  the  plume  is  embedded.  If 


Table  1 

Description  of  the  main  physical  features  of  the  initial  and  final  stages  of  Deep  Convection 
Q  NL  LA  PM 


TM 


Initial  stages  No 
Final  stages  Yes 


Yes  Yes  No 

No  No  Yes 


Yes  Yes 

No  Yes 


Non-locality  (NL),  lateral  advection  due  to  baroclinic  instability  (LA),  plume  model  (PM)  and  turbulence  model  (TM). 
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the  downward  plume’s  velocity  is  wd  and  the  updraft’s  velocity  is  wu,  mass  conservation  demands 
that: 

(1  -  a) /a  =  wu/wd  (2d) 

Since  in  the  initial  stages  of  development  downward  plumes  have  velocities  some  ten  times  larger 
than  the  updrafts,  Eq.  (2d)  implies  again  (2c).  Thus,  both  assumptions  of  the  OPPS  model  rep¬ 
resent  the  same  physical  picture  of  a  plume  in  the  initial  stages  of  development  whereas  a  tur¬ 
bulence  model  can  describe  all  stages  provided  it  is  non-local,  that  is,  capable  of  incorporating  the 
velocity  field  skewness  Sw. 

(4)  Deep  convection  is  not  a  purely  vertical  process.  The  initial  and  final  stages  are  characterized 
not  only  by  different  values  of  the  skewness  Sw  but  also  by  an  additional  important  feature  (Jones 
and  Marshall,  1993;  Visbeck  et  al.,  1996;  Marshall  and  Schott,  1999;  Herbaut  and  Marshall,  2002; 
Lavender  et  al.,  2002).  In  the  initial  phases,  the  mixed  layer  deepens  at  a  rate  in  accord  with  that 
given  by  vertical  mixing  models  (Visbeck  et  al.,  1996).  However,  when  the  plume  has  deepened  to 
a  distance  £  such  that  £  >  ^(rot),  rotation  takes  over  and  a  baroclinic  instability  sets  in.  The 
deepening  process  is  arrested  since  baroclinic  eddies  transport  fluid  masses  away  laterally.  Thus, 
in  the  later  stages  of  DC,  lateral  advection  dominates  over  vertical  transport  (see  Fig.  6  of  Visbeck 
et  al.,  1996).  Since  Lavender  et  al.  (2002,  Fig.  13)  have  shown  that  the  heat  flux  carried  by  the 
vertical  plumes  cannot  balance  the  winter  surface  heat  loss,  a  considerable  fraction  of  the  heat  flux 
must  be  carried  by  “features”  that  have  a  life  time  longer  than  that  of  the  plumes  0(hours).  In 
Table  1  we  highlight  the  relative  importance  of  the  different  processes  in  the  initial  and  final  stages 
of  deep  convection. 

Since  a  plume  model  is  applicable  mostly  in  the  initial  stages  of  deep  convection,  it  may  be  difficult 
to  connect  it  with  the  initiation  of  baroclinic  instabilities  that  occur  mostly  in  the  final  stages  (LA 
and  PM  columns).  On  the  other  hand,  a  turbulence  model  valid  for  arbitrary  values  of  skewness  S„ 
and  rotation,  is  applicable  to  all  stages  of  deep  convection  (TM  column).  In  conclusion,  a  reliable 
model for  deep  convection  requires  the  combination  of  a  vertical  mixing  turbulence  model  with  rotation 
and  non-locality,  together  with  a  mesoscale  parameterization  to  describe  lateral  advection. 


2.  Models  of  deep  convection 

2.1.  Convective  adjustment 

Convective  Adjustment  (CA)  was  proposed  by  Bryan  (1969)  as  a  temporary  model  but  it  is  still 
in  use  today.  It  has  been  recently  analyzed  by  Marotzke  (1991)  and  Marshall  and  Schott  (1999).  In 
the  latter  study,  it  was  stressed  that  CA  requires  an  instantaneous  adjustment  of  the  water  column 
while  in  reality  the  time  of  mixing  tmix  is  finite  (of  the  order  of  12  h  or  so).  Stated  differently,  the 
CA  adjustment  scheme  corresponds  to  an  extremely  large  vertical  diffusivity  Kv  whereas  one 
expects  that: 

K\  &  h  / tm jx  ~  hWp  (3) 

Using  wp  —  5  cm/s  and  h  =  1  km,  one  obtains  Kv  ~  50  m2s_1  (Send  and  Marshall,  1995).  Mar¬ 
otzke  (1991)  used  Kv  =  1  m2s_1  while  Klinger  et  al.  (1996)  used  Kv  =  (10-50)  m2s_1.  From  the 
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viewpoint  of  the  consequences  of  the  CA  scheme,  Bryan  (1986)  pointed  out  that  the  standard 
GFDL  scheme  leads  to  the  collapse  of  the  meridional  circulation.  On  the  other  hand  the  collapse 
is  avoided  if  one  employs  a  complete  mixing  scheme  and/or  a  treatment  of  DC  as  a  vertical 
diffusion  (Marotzke,  1991).  Reviewing  the  problem,  Paluszkiewicz  and  Romea  (1997)  concluded 
that  the  CA  schemes  do  not  produce  realistic  vertical  density  structure,  do  not  create  the  correct 
quantity  of  deep  water,  and  do  not  use  a  time  scale  of  adjustment  in  agreement  with  tracer  ages  or 
observations.  Rim  and  Stossel  (2001)  have  pointed  out  that  CA  yieldsexcessive^ 
upwells  too  much  heat  which  in  turn  yields  too  little  sea-ice  coven  Finally,  Yin  and  Sarachik 
119941  suggested  and  tested  a  scheme  that  treats  DC  as  a  vertical  diffusion. 

It  would  be  reassuring  if  a  turbulence  model  could  reproduce  the  values  of  the  vertical  diffu- 
sivities  used  in  these  heuristic  models.  As  we  shall  show  m  Fig.  4,  this  is  indeed  the  case. 


2.2.  Turbulence  models 

Conspicuously  absent  in  the  treatment  of  deep  convection  have  been  turbulence  models,  as 
pointedout  by  several  authors.  Killworth  (1989)  has  stated  that  “while  models  for  the  mixed  layer 
dvnamic  are  legion  those  for  deep  convection  are  rare”;  Maxworthy  (1997)  wrote  that:  thus  far, 
turbulence  schemes  have  tended  to  be  restricted  to  surface  andi  local  area  modeling’  ^  that  one 
wants  to  hope  that  more  complex  closure  schemes  will  be  tried  on  problems  like  deep  convection 
while  Marshall  and  Schott  (1999)  state  that:  “an  effort  should  be  made  to  apply  them  to  these 
deep  convective  flows”  and  that  “the  challenge  for  the  future  is  to  transform  the  insights  about 

deep  convection  into  a  parametric  representation”. 

This  paper  is  an  attempt  to  fill  this  gap.  Specifically,  we  assess  the  MY  KPP  and  GISS  models 
for  the  vertical  mixing.  Lateral  advection  is  modeled  with  the  parameterization  of  Gent  and 
McWilliams  (1990).  The  data  on  deep  convection  are  from  the  Labrador  Sea  (Lavender  e  ., 
^  In  Us  present  formulation,  the  GISS  model  is  local  (valid  in  the  final  stages  of  evolu  ion 
third  column  in  Table  1)  and  has  no  rotation.  Thus,  we  expect  that  the  model  will  overestimate  the 
convective  depth  as  the  results  presented  below  will  indeed  show.  In  a  second  study,  we  shall  use 
the  GISS  loca? model  (plus  the  GM  model)  but  with  the  inclusion  of  rotation  m  the  vertical  mixing 
model.  In  a  third  study,  we  shall  include  non-locality  in  the  GISS  model  so  as  to  account  for  both 
initial  and  final  stages  of  deep  convective  processes  (lateral  advection  is  still  treated  with  the  GM 
^el)  In  a  fourth  study,  we  shall  use  the  GISS  model  employed  in  the  third  study,  but  the 
meiscale  model  will  no  longer  be  the  Gent-McWilliams  model  but  a  new  one  (Canuto  and 
Dubovikov,  submitted  for  publication)  that  has  been  developed  and  is  presently  being  tested 
against  data  from  an  eddy  resolving  ocean  code. 


3.  The  GISS  turbulence  model:  general  features 

This  vertical  mixing  model  was  constructed  using  the  Reynolds  stress  method  (RSM,  Canuto 
etll  2001a  2002;  did  as  Cl, 2;  Cheng  et  al.,  2002).  It  includes  a  velocity  field  and  two  active 
scalar  fields,  temperature  and  salinity.  The  model  can  be  time  dependentorstationaryoca  or 
non-local,  with  or  without  rotation.  In  the  stationary  and  local  case,  the  GISS  model  becomes  an 
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algebraic  set  of  equations  whose  solution  is  analytical.  The  resulting  vertical  turbulent  fluxes  of 
momentum,  heat  and  salt  are  given  by 


W 
9 z 


mV  =  ~Km  — ,  wT  =  -Kh 


0 T 
dz  ’ 


w>s'  =  -Ks  — 


es 

dz 


(4a) 


Here,  U,  T  and  S  are  the  mean  velocity,  temperature  and  salinity.  The  vertical  eddy  diffusivities  Ka 
(a  stands  for  momentum,  heat  and  salinity)  can  be  written  in  two  representations 


Kai  =  rale2z,  Ka 2  =  ra2sN~2 


(4b) 


where  the  mixing  efficiencies  ra  are  given  by 

ral  =  cS^Zf1 ,  2r<a  =  S«(Nx)2  (4c) 

Here,  i  is  the  mixing  length,  e  is  the  dissipation  of  turbulent  kinetic  energy,  t  is  the  eddy  turnover 
time,  N  is  the  Brunt-Vaisala  frequency,  Z  is  the  mean  shear  and  c  is  a  constant.  The  structure 
functions  Sa  are  analytic  functions  of  the  Richardson  number  Ri  =  N2/Z2  and  of  the  density  ratio 
Rp  =  pdS/dz(adT/dz)~l  where  a  and  /?  are  the  thermal  expansion  and  haline  contraction  coeffi¬ 
cients.  Once  the  heat  and  salt  diffusivities  are  known,  the  passive  scalar  Kc  and  mass  diffusivities 
Kp  are  also  known  since 

Kc  =  (Kh  +  RpKs)(l  +Rp)~\  Kp  =  {Kh-RPKS){  1  -  Rp)~l  (4d) 

The  physical  interpretation  of  the  model  is  as  follows.  Ka i  is  the  contribution  due  to  wind  shear 
while  Ka  2  is  due  to  internal  wave  breaking  if  no  convection  is  present.  Since  these  processes  coexist, 
the  diffusivity  Ka  is  the  sum  of  the  two  contributions 

Ka=Kal+Ka 2  (4e) 

In  the  ML,  where  wind  shear  is  much  larger  than  the  internal  wave  breaking  contribution,  the  first 
term  in  (4e)  dominates;  below  the  ML,  where  wind  shear  becomes  inefficient  (large  Ri),  the  second 
term  in  (4e)  dominates  and  thus  a  smooth  transition  of  processes  is  assured.  When  N2  <0 
(convective  regime),  the  contribution  to  Ka2  by  internal  wave  breaking  becomes  negligible  and  Ka2 
is  predominantly  due  to  convection. 


4.  Previous  tests  of  the  GISS  model 

Mixed  layer  and  thermo  dine.  The  first  representation  in  (4b)  was  used  in  Cl, 2  and  by  Burchard 
and  Bolding  (2001)  and  Burchard  and  Deleersnijder  (2001).  Below  the  mixed  layer  and  in  the 
absence  of  deep  convection,  the  most  likely  source  of  mixing  is  the  breaking  of  internal  waves. 
Using  the  second  representation  (4b),  the  value  of  eN~2  (N2  >  0)  due  to  internal  waves  breaking 
was  taken  from  (Polzin,  1996). 

Convective  efficiencies  ra.  Since  their  original  introduction  (Osborn  and  Cox,  1972),  many 
discussions  have  concentrated  over  the  values  of  these  convective  efficiencies  which,  without  a 
turbulence  model,  remain  adjustable  parameters.  Several  questions  have  remained  without  an 
answer,  for  example:  (1)  are  the  ra2  for  momentum,  heat  and  salt  the  same?  (2)  is  it  correct  to 
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assume  Ta2  =  1/4  as  often  done,  when  large  eddy  simulation  (LES)  data  (Wang  et  al.,  1996)  yield 
the  much  larger  value  Ta2  ~  1  ? 

The  GISS  model  computes  the  TVs  which  were  shown  to  increase  strongly  near  Rp  =  0.6,  a 
prediction  that  is  consistent  with  the  recent  observations  of  a  much  larger  mixing  at  Barbados 
( Rp  =  0.6)  than  at  the  NATRE  location  (Ledwell  et  al.,  1998)  where  Rp  =  0.56.  The  GISS  model 
also  predicts  values  of  Ta 2  up  to  unity,  in  agreement  with  LES  data  (Wang  et  al.,  1996)  thus 
casting  doubts  about  the  universality  of  the  commonly  used  values  0.2-0.25. 

NATRE  data.  KSih  have  often  been  assumed  to  be  the  same  but  the  NATRE  measurements  have 
shown  that  Kh  ^  Ks.  Several  authors  have  investigated  the  influence  of  double  diffusive  processes 
in  OGCMs  (Zhang  et  al.,  1998;  Merryfield  et  al.,  1999).  Some  authors  treated  the  ratio  Kh/Ks  as  a 
free  parameter,  others  (Large  et  al.,  1994)  employed  the  relation  Ks  =  1 .43 (however,  the 
measurements  show  that  Ks  =  l.43KhR~l).  In  general,  these  studies  conclude  that  double  diffusion 
alters  substantially  the  regional  distribution  of  T  and  S  and  that  salt  fingers  are  more  widespread 
than  diffusive  convection.  Compared  with  the  Kh  =  Ks  case,  deep  temperature  and  salinity  become 
larger,  the  polar  heat  transport  is  lowered  and  the  meridional  overturning  is  reduced.  Given  the 
heuristic  nature  of  the  relation  Ks/Kh  employed  thus  far,  it  is  difficult  to  compare  and/or  assess  the 
reliability  of  the  ensuing  results. 

Tests  without  an  OGCM.  Since  at  400  m  depth  the  diffusivities  are  only  functions  of  Rp,  the 
values  of  Rp(z )  measured  at  the  NATRE  site  (Ledwell  et  al.,  1998)  were  used  to  compute  the 
concentration  and  mass  diffusivities  from  Eq.  (4d).  Comparison  of  the  GISS  model  results  with 
the  NATRE  data  was  presented  in  Fig.  9  of  C2.  Unpublished  results  from  data  at  Barbados 
(R.W.  Schmitt,  Hawaii  Ocean  Sciences  Meeting,  Feb.  2002)  of  a  much  larger  mixing  at  Barbados 
than  at  NATRE,  are  also  in  qualitative  agreement  with  the  predictions  of  the  model. 

Tests  with  an  OGCM.  The  GISS  model  was  used  in  a  coarse  resolution  GFDL-type  ocean  code 
and  the  heat,  salt,  mass  and  concentration  diffusivities  at  the  NATRE  location  were  computed 
and  shown  to  reproduce  the  NATRE  data  quite  closely.  The  775  profiles  in  different  ocean’s  sites 
were  also  computed  and  in  general  they  reproduce  the  Levitus  data  well.  Particularly  significant  is 
the  better  agreement  of  the  T  vs.  z  profile  in  the  Arctic  Ocean. 


5.  Modeling  deep  convection 

5.1.  3°  x  3°  resolution 

We  have  employed  the  NCAR-CSM  global  ocean  model  (Large  et  al.,  1997;  for  details,  see 
Appendix  A)  with  a  3°  x  3°  resolution.  The  model  requires  a  vertical  mixing  scheme  and  employs  a 
GM  model  to  parameterize  the  mesoscale  eddies.  We  ran  two  global  simulations  with  two  dif¬ 
ferent  vertical  mixing  schemes,  the  KPP  model  (Large  et  al.,  1994)  and  the  GISS  model  discussed 
above.  The  monthly  surface  forcing  calculation  and  all  other  aspects  of  the  two  OGCM  runs  were 
identical.  To  compare  the  model  results  with  observations,  in  Fig.  1  we  reproduce  the  results  of 
Lavender  et  al.  (2002).  Measurements  show  that  the  deepest  convection,  DC,  one  with  mixed- 
layers  deeper  than  0.8  km,  is  confined  to  the  fairly  small  Lavender,  Davis  and  Owens  (LDO) 
region  in  the  western  Labrador  Sea,  Fig.  1.  A  similar  result  is  presented  in  Fig.  12  of  Pickart  et  al. 
(2002).  In  Figs.  2  and  3,  we  exhibit  the  GISS  and  KPP  simulation  points  found  to  yield  deepest 
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Fig.  1.  Mixed  layer  depths  as  from  Lavender  et  al.  (2002). 


Fig.  2.  DC  (in  km)  as  from  the  GISS  mixing  model. 

convection  (the  results  are  snapshots  taken  at  March  1  of  the  127th  year).  The  GISS  mixing  model 
reproduces  the  deepest  convection  in  the  LDO  region  while  the  KPP  model  does  not.  Within  the 
LDO  region,  the  KPP  mixed  layers  are  not  deeper  than  20  m.  Deep  convection  in  the  GISS  model 
is  evident  in  both  the  flat  vertical  temperature  profile  and  large  heat  diffusivity  down  to  1.2  km 
(Fig.  4).  The  GISS  model  heat  diffusivity  in  this  column  increases  from  0.04  m2s_1  near  the 
surface  to  a  maximum  of  6  m2  s-1  at  400  m.  In  other  columns,  we  find  maximum  diffusivities  up  to 
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Fig.  3.  DC  (in  km)  as  from  the  KPP  mixing  model. 


Fig.  4.  Profiles  of  the  temperature  and  heat  diffusivity  (inset)  at  52.2W,  57. 5N  within  the  box  of  Fig.  1. 


10  m2s-1.  The  maximum  depth  attained  in  the  GISS  model  falls  within  the  observational  range 
while  that  for  the  KPP  model,  which  occurs  outside  the  LDO  box,  is  larger  than  observed.  The 
temperature  and  diffusivity  profiles  for  the  model  point  52.2W,  57. 5N  nearest  the  center  of  the 
“deep  convection  box”,  are  shown  in  Fig.  4.  In  the  LDO  region,  the  GISS  model  diffusivity  lies 
between  1  m2  s“‘  used  by  Marotzke  (1991)  and  10-50  m2  s-1  used  by  Klinger  et  al.  (1996).  At  the 
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points  where  the  KPP  model  produces  deep  convection,  the  diffusivities  are  of  the  same  order  of 
magnitude  as  the  GISS  model,  though  generally  somewhat  larger.  We  also  present  globally  and 
annually  averaged  temperature  and  salinity  profiles  vs.  Levitus  data  from  the  126th  model  year. 
Figs.  5  and  6  ought  to  be  compared  with  Figs.  12  and  13  of  C2.  The  KPP  model  results  are  almost 
unchanged  (although  very  slightly  fresher)  while  the  GISS  model  results  are  significantly  im¬ 
proved.  The  cool  bias  of  the  GISS  model  vs.  observations  at  and  below  1  km  has  been  almost 
eliminated  and  the  salinity  bias  between  1  and  3  km  has  been  considerably  reduced. 


5.2.  Sensitivity  tests 

To  investigate  the  sensitivity  of  the  location  of  deep  convection  to  model  configuration,  in 
particular  the  landmask,  we  ran  two  tests  with  geographies  different  from,  though  less  realistic 
than,  the  standard  case  used  for  the  results  described  above.  The  addition  and  subtraction  of  land 
points  was  used  as  a  simple  way  to  perturb  the  ocean  model  in  order  to  assess  the  extent  to  which 
the  comparison  between  the  different  vertical  mixing  schemes  would  be  affected  by  an  extraneous 
factor.  In  each  test  we  altered  two  points  of  the  landmask  while  keeping  everything  else  the  same. 

In  test  A  (“+LAND”)  we  added  land  to  northern  Labrador.  The  two  model  points  at  (59.4W, 
57.5N)  and  (59.4W,  59.38N)  were  converted  from  16  level  (1.6  km  deep)  ocean  points  into  land 
points.  The  results  of  test  A  are  shown  in  Fig.  7a  and  b.  The  GISS  model  has  deepest  convection 


Fig.  5.  Global  annually  averaged  temperature  vs.  depth  for  both  mixing  models  compared  to  the  Levitus  data.  As¬ 
terisks,  GISS  model;  diamonds,  KPP  model;  solid  line,  Levitus. 
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Fig.  6.  Global  annually  averaged  salinity  vs.  depth  for  both  mixing  models  compared  to  the  Levitus  data.  Asterisks, 
GISS  model;  diamonds,  KPP  model;  solid  line,  Levitus. 


Fig.  7.  (a)  DC  (in  km)  as  from  the  GISS  mixing  model.  Case  +  LAND  (land  added  to  northern  Labrador),  (b)  DC 
(in  km)  as  from  the  KPP  mixing  model.  Case  +  LAND  (land  added  to  northern  Labrador). 


in  fewer  spots,  one  of  these  still  being  the  southern  point  of  the  box  where  the  observations  in 
Lavender  et  al.  (2002)  exhibit  it,  but  unlike  the  standard  case,  there  is  no  deepest  convection  in  the 
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Fig.  8.  (a)  DC  (in  km)  as  from  the  GISS  mixing  model.  Case  +  SEA  (sea  extended  southward  into  Labrador),  (b)  DC 
(in  km)  as  from  the  KPP  mixing  model.  Case  +  SEA  (sea  extended  southward  into  Labrador). 

north  of  the  box.  The  KPP  also  has  fewer  convection  spots,  but  they  still  extend  over  a  wide 
region  and  yet,  no  deep  convection  appears  within  the  box. 

In  test  B  (“+SEA”)  we  extended  the  sea  southward  into  Labrador.  The  two  model  points  at 
(59.4W,  55.51N)  and  (59.4W,  53.42N)  were  converted  from  land  points  into  16  level  ocean  points. 
The  results  of  test  B  are  shown  in  Fig.  8a  and  b.  The  GISS  model’s  deepest  convection  points  are 
now  slightly  deeper  on  average  and  it  now  has  deepest  convection  in  the  north  of  the  box,  but  not 
in  the  south.  The  KPP  model’s  deepest  convection  moves  somewhat  to  the  northeast  but  there  is 
still  no  deep  convection  inside  the  LDO  box.  In  both  tests  A-B,  no  KPP  mixed  layer  inside  the 
LDO  box  was  deeper  than  20  m.  Although  these  results  show  that  addition  or  subtraction  of 
landpoints  affects  DC,  the  global  results  have  remained  largely  unchanged.  The  ocean  basin 
profiles  for  the  KPP  and  GISS  models  in  the  test  A  and  B  cases  are  very  similar  to  those  obtained 
with  the  standard  landmask  and  for  this  reason  we  do  not  reproduce  them  here. 

While  adding  or  taking  away  land  degraded  the  GISS  model’s  performance  vis-a-vis  the  more 
realistic  standard  landmask  case,  it  did  not  improve  the  KPP  model’s  performance.  Furthermore, 
in  both  cases  the  GISS  model  remained  more  realistic  in  having  some  deepest  convection  inside 
the  observed  box.  The  qualitative  conclusion,  robust  under  change  of  the  landmask  for  the 
NCAR  CSM  ocean  model  tests,  is  that  both  models  produce  deepest  convection  in  the  Labrador 
Sea  area  over  a  larger  region  than  observed  but  that  only  the  GISS  model’s  deepest  convection 
overlaps  the  observations. 

5.3.  1°  x  1°  resolution 

In  this  case,  the  vertical  mixing  scheme  is  the  MY2.5  model  (S.  Hakkinen,  private  communi¬ 
cation)  and  it  was  used  in  an  ocean  model  of  the  northern  hemisphere  (cr-coordinates  with  20 
levels  with  higher  resolution  near  the  surface,  with  monthly  forcing,  Hakkinen,  1999,  2001).  In 
Fig.  9a  and  b  we  present  the  model  results.  Fig.  9a  corresponds  to  a  forcing  from  climatology 
while  in  Fig.  9b  the  forcing  corresponds  to  the  year  1998.  The  model  produces  deepest  convection 
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Fig.  9.  (a)  DC  (in  km)  as  from  the  MY2.5  mixing  model.  Climatological  forcing,  (b).  DC  (in  km)  as  from  the  MY2.5 
mixing  model.  The  forcing  corresponds  to  1998. 

within  the  LDO  region  with  depths  similar  to  those  observed  but  such  regions  extend  west  and 
much  further  north  and  south  than  observed.  The  depths  in  Fig.  9b  are  somewhat  shallower  than 
those  in  Fig.  9a  but  the  extent  of  the  deepest  convection  region  is  almost  as  large  as  before. 

5.4.  113°  x  1/3°  and  1/12°  x  1/12°  resolutions 

Two  more  simulations  were  performed  with  the  hybrid  coordinate  ocean  model  (HYCOM; 
Bleck,  2002).  HYCOM  is  isopycnal  in  the  open,  stratified  ocean,  but  using  the  layered  continuity 
equation  makes  a  dynamically  smooth  transition  to  a  terrain-following  (sigma)  coordinate  in 
shallow  coastal  regions,  and  to  z-level  coordinates  in  the  mixed  layer  and/or  unstratified  seas.  In 
the  North  Atlantic  HYCOM,  the  vertical  coordinates  are  configured  to  allow  high  vertical  res¬ 
olution  (order  15  z-level  coordinates  in  the  top  200  m  during  the  winter)  in  the  mixed  layer, 
thereby  allowing  the  use  of  the  KPP  closure  scheme.  The  model  includes  the  entire  North  Atlantic 
(70°N  to  28°S),  but  only  the  Labrador  Sea  is  discussed  here.  There  are  26  layers  in  the  vertical, 
and  horizontal  resolution  is  1/3°  (discussed  here),  and  1/12°  (discussed  below).  The  MODAS 
climatology  (Fox  et  al.,  2002)  was  used  to  initialize  the  model  and  to  parameterize  the  meridional 
overturning  circulation  via  relaxation  to  temperature  and  salinity  in  3°  buffer  zones  at  the 
northern  and  southern  boundaries.  A  monthly  climatology  formed  from  the  European  Centre  for 
Medium  Range  Forecasting  (ECMWF)  10  m  reanalysis  is  used  for  surface  forcing.  Bulk  formulae 
were  used  to  calculate  the  surface  heat  flux  from  the  net  radiation  flux,  the  air  temperature,  the 
specific  air  humidity,  the  scalar  wind  speed,  and  the  model’s  surface  (layer  1)  temperature.  Using 
the  model’s  SST  (instead  of  climatological  SST)  tends  to  self  correct  the  net  surface  heat  flux, 
i.e.,  if  the  SST  is  too  warm  the  heat  flux  will  decrease  and  vice  versa. 

For  ease  of  comparison,  the  results  corresponding  to  March  1  from  the  HYCOM  simulations 
were  binned  in  2.5°  squares  centered  on  the  same  grid  points  as  the  NCAR-CSM  results  shown  in 
Figs.  2  and  3.  The  maximum  mixed  layer  depth  within  each  bin  that  exceeds  800  m  is  shown  in 
Fig.  10a,  and  the  mean  mixed  layer  depth  of  each  bin  (greater  than  800  m)  is  shown  in  Fig.  10b. 
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Fig.  10.  (a)  The  maximum  (larger  than  0.8  km)  mixed  layer  depths  (in  km)  (within  a  3°  box)  predicted  by  KPP 
(HYCOM  1/3°  x  1/3°  resolution  ocean  code).  Forcing  is  with  monthly  ECMWF  winds  and  thermal  fields,  (b)  Same  as  in 
(a)  but  for  the  mean  mixed  layer  depths  (within  a  3°  box). 

In  HYCOM,  the  mixed  layer  depth  is  defined  as  the  minimum  depth  where  the  aQ  change  from 
the  surface  downward  corresponds  to  a  0.2  °C  temperature  change.  Overall,  the  maximum  mixed 
layer  depths  are  deeper  than  those  observed  as  shown  in  Fig.  1,  and  the  distribution  of  the 
maxima  is  not  consistent  with  the  observations.  There  are  only  two  mean  mixed  layer  depth  bins 
^  800  m  (deepest  convection),  but  their  locations  are  consistent  with  the  outlined  LDO  region 
(Fig.  10b).  In  HYCOM  with  KPP,  the  deep  mixed  layers  are  often  spatially  uncorrelated,  i.e., 
there  can  be  a  1 500  m  deep  mixed  layer  and  100  m  deep  mixed  layer  at  adjacent  grid  points.  This  is 
apparently  caused  by  the  magnification  of  small  scale  horizontal  numerical  noise  by  the  mixed 
layer  model.  This  may  be  due  to  the  one-dimensional  nature  of  KPP,  but  other  factors,  such  as 
how  mesoscale  ocean  features  respond  to  large  scale  atmospheric  forcing,  may  also  play  a  role. 
The  result,  however,  is  that  on  this  particular  (climatological)  day,  the  mean  mixed  layer  depths 
are  heavily  influenced  by  a  larger  number  of  “shallow”  mixed  layer  depths  within  each  bin, 
thereby  decreasing  the  number  of  bins  that  meet  the  800  m  minimum  threshold  (the  region  of 
deepest  convection  as  defined  earlier). 

5.5.  1/12°  x  1/12° 

The  1/3°  HYCOM  model  described  in  the  previous  section  was  run  with  an  horizontal  reso¬ 
lution  of  1/12°  (about  4.7  km  at  58°N).  The  maximum  and  mean  mixed  layer  depths  in  the  2.5° 
bins  are  shown  in  Fig.  1  la  and  b,  respectively.  As  with  the  1/3°  results,  there  are  more  maximum 
mixed  layer  depths  than  mean  mixed  layer  depths,  for  the  reasons  described  above.  However,  in 
the  1/12°  model,  the  mixed  layer  depths  are  much  deeper  than  in  the  1/3°  case.  The  range  of  the 
mean  mixed  layer  depths  (Fig.  lib)  is  0.8-1. 7  km.  These  depths  are  consistent  with  the  obser¬ 
vations,  although  most  of  the  deep  convection  sites  lie  to  the  south-east  of  the  observed  deepest 
convection  sites  (Fig.  1).  The  maximum  mixed  layer  depths  (Fig.  11a)  are  much  deeper  than  those 
observed,  with  several  of  the  order  of  3  km.  They  are  also  more  widely  distributed  than  observed. 


90 


V.  M.  Canuto  et  at  /  Ocean  Modelling  7  ( 2004)  75-95 


Fig.  11.  (a)  Same  as  in  Fig.  10a  with  a  1/12°  x  1/12°  resolution,  (b)  Same  as  in  Fig.  10b  with  a  1/12°  x  1/12°  resolution. 

This  may  be  due  to  limitations  in  how  KPP  parameterizes  deep  convective  turbulent  boundary 
layers  and/or  due  to  other  factors  such  as  HYCOM’s  hybrid  coordinate  or  more  inertial  currents 
that  give  rise  to  higher  current  shears. 


6.  Improvements  on  the  GISS  model 

6.1.  Rotation 

As  we  have  already  discussed,  this  is  the  first  in  a  series  of  studies  of  DC  that  we  have  planned. 
Here,  we  have  assessed  the  performance  of  a  local  model  for  vertical  mixing  plus  a  GM  model  for 
the  mesoscale  eddies.  The  inclusion  of  rotation  in  the  vertical  mixing  model  will  be  as  follows.  In 
the  one-point  closure  model  that  we  employ,  rotation  enters  in  two  ways:  through  the  Coriolis 
force  that  comes  in  linearly  in  the  equations  and  into  the  pressure-correlations  (Canuto  et  al., 
2001a).  Such  pressure  correlations  depend  on  shear  2 Sy  =  Utj  +  Ujj  and  vorticity  2  Vy  =  Uij  —  Ujj, 
where  U,  is  the  mean  velocity  field  and  Uu  =  dUi/dxj.  In  the  presence  of  Q,  the  Vv  change  to 

Vij  -  V*  =  V,j  -  eiJkQk  M 

where  eijk  is  the  totally  antisymmetric  tensor.  Since  Vv  enters  in  all  turbulence  equations,  both  the 
heat  and  the  momentum  fluxes  (Reynolds  stresses)  are  affected  by  rotation  Q  =  (0,  Qy,  Qz),  where 
Qyf  =  Q (cos  cp,  sirup),  cp  being  the  latitude.  The  dynamic  equations  for  the  heat  fluxes  and  the 
Reynolds  stresses  in  the  presence  of  £2  have  recently  been  derived.  In  the  local  limit,  they  become 
an  algebraic  system  that  we  have  already  solved.  As  an  example,  the  first  of  (4a)  now  contains 
additional  terms  of  the  form: 

AQ cos  cp  +  BQ  sin  cp 

where  cp  is  the  latitude.  A  similar  structure  also  enters  the  heat  flux. 


(5b) 
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6.2.  Non-locality 

For  a  vertical  mixing  model  to  be  valid  both  at  the  beginning  and  at  the  end  of  the  DC  process, 
one  must  allow  arbitrary  values  of  the  skewness  Sw  which,  as  Eq.  (2a)  shows,  is  a  third-order 
moment  (TOM).  Non-local  terms  appear  in  both  the  equation  for  the  heat  flux  as  well  as  in  that  of 
the  Reynolds  stresses.  Before  we  discuss  how  to  model  the  TOMs,  we  note  that  in  the  non-local 
case,  Eq.  (4a)  changes  to 

uW  =  —KmdU/dz  +  ym,  wT  =  -KhdT/dz  +  yh  (6a) 

While  these  expressions  are  not  new,  previous  authors  (e.g.,  Large  et  al.,  1994)  had  to  employ 
heuristic  models  for  the  y’s  and  for  lack  of  data,  had  to  take  ym  =  0.  On  the  other  hand,  a  tur¬ 
bulence  model  provides  the  y’s  in  terms  of  the  TOMs,  e.g.: 

ft  ft 

y  h~—W2r+A—  WT1  (6b) 

dz  oz 

Physically,  the  first  term  represents  the  flux  of  the  heat  flux  while  the  second  term  represents  the 
flux  of  the  temperature  variance.  The  presence  of  the  y’s  shows  that  turbulence  generates  and 
transports  fluxes.  Thus,  a  flux  at  a  given  site  need  not  have  been  generated  locally;  it  may  be  the 
result  of  the  transport  of  a  flux  generated  elsewhere,  that  being  the  meaning  of  non-locality. 
Heuristic  expressions  for  the  TOMs  were  employed  for  many  years  until  Moeng  and  Wyngaard 
(1989)  showed  that  they  provide  a  very  poor  representation  of  the  TOMs  obtained  from  LES. 
Motivated  by  that  result,  a  new  expression  for  the  TOMs  was  derived  (Canuto  et  al.,  1994)  by 
solving  their  dynamic  equations.  The  results  were  tested  with  good  results  against  three  sets  of 
LES  data  for  a  convective  PBL.  Recently,  a  new  expression  for  the  TOMs  has  been  derived 
(Canuto  et  al.,  2001b)  which  is  more  compact  than  the  1994  one.  When  compared  with  a  new  set 
of  LES  data  for  the  planetary  boundary  layer  PBL  (Mironov  et  al.,  2000),  it  yields  results  that  are 
superior  to  those  in  1994.  However,  ocean  TOMs  are  not  the  same  as  those  for  the  PBL.  In  fact, 
since  earth’s  rotation  influences  the  dynamics  of  ocean  DC,  the  TOMs  must  depend  on  Q.  To 
derive  the  new  TOMs(fi),  we  shall  employ  the  methodology  used  in  the  two  previous  models  for 
the  TOMs.  The  new  TOMs(fi)  will  first  to  be  tested  against  new  LES  data  that  include  rotation. 

6.3.  Mesoscale  modeling 

In  principle,  mesoscale  eddies  affect  both  the  heat/salt  fluxes  as  well  as  momentum  fluxes 
(Reynolds  stresses).  Since  the  Gent  and  McWilliams  model  (1990)  was  devised  for  scalar  fluxes, 
the  Reynolds  stresses  are  still  modeled  through  a  Fickian  diffusion.  Improvements  of  the  GM 
model  are  difficult  since  the  model  was  not  derived  from  a  set  of  dynamic  equations.  However, 
starting  from  the  Navier-Stokes  equations  and  the  equation  for  a  passive  scalar,  Canuto  and 
Dubovikov  (submitted  for  publication)  have  recently  derived  a  model  whereby  the  mesoscale 
effect  on  both  the  mean  momentum  equations  and  on  the  equations  for  scalars  can  be  computed. 
Some  key  new  features  of  the  new  model  are  as  follows: 

(1)  the  bolus  velocity  contains  two  terms,  the  first  of  which  is  the  downgradient  of  the  potential 
vorticity  rather  than  the  mean  thickness  as  suggested  by  several  authors  (e.g.,  Treguier,  1999), 


92 


V.M.  Canuto  et  al.  /  Ocean  Modelling  7  (2004)  75-95 


(2)  the  second  term  depends  on  the  mean  velocity  field  u  and  is  not  present  in  any  of  the  param- 
eterizations  discussed  thus  far  in  the  literature.  The  existence  of  a  new  term  was  first  suggested 
by  Bryan  et  al.  (1999)  but  no  specific  form  was  derived  for  it. 

(3)  since  u  varies  with  depth,  the  new  term  causes  a  shearing  force  that  erodes  the  coherence  of  the 

eddy  field,  . 

(4)  thus  far,  the  mesoscale  diffusivity  Km  was  treated  either  as  an  adjustable  parameter  or  dealt 

with  using  a  heuristic  model  (Visbeck  et  al.,  1997).  The  new  model  expresses  Km  in  terms  of 
the  large  scale  (resolved)  fields  and  thus  /cm  can  be  computed.  It  is  predicted  to  be  larger  at 
the  surface  than  in  the  interior,  in  agreement  with  recent  heuristic  models  (Visbeck  et  al., 
1997;  Karsten  and  Marshall,  2002). 


7.  Conclusions 

Since  our  goal  was  to  investigate  models  for  deep  convection  to  be  used  in  coarse  resolution 
OGCMs,  we  have  analyzed  three  vertical  mixing  models.  Considering  them  in  the  historical  order 
in  which  they  appeared,  we  have  discussed:  (a)  MY2.5,  (b)  KPP  and  (c)  GISS  model. 

(a)  as  discussed  in  Canuto  et  al.  (2001a,b,  2002),  in  stably  stratified  situations,  the  MY  models 
produce  too  shallow  mixed  layers,  a  point  first  made  by  Martin  (1985).  Technically,  this  is 
due  to  the  fact  that  the  MY  models  predict  a  critical  Richardson  number  Ri(cr)  =  0.19  which 
is  five  times  smaller  than  the  value  needed  to  obtain  correct  mixed  layer  depths.  In  an  unstably 
stratified  configuration,  such  as  the  one  characterizing  deep  convection,  the  MY  models  per¬ 
form  better,  as  shown  in  Fig.  9.  However,  since  a  mixing  model  is  supposed  to  be  applicable 
everywhere,  the  shortcomings  of  the  MY  models  in  stable  situations  must  still  be  amended. 

(b)  the  KPP  model  was  constructed  primarily  to  alleviate  the  “too  little  mixing”  problem  of  the 
MY  models.  It  is  interesting  to  compare  its  performance  vis-a-vis  the  MY  models:  KPP  yields 
correct  mixed  layer  depths  in  stable  situations  but,  at  coarse  resolution,  it  produces  too  little 
mixing  in  the  LDO  region  (Figs.  3  and  4).  However,  as  the  resolution  increases,  the  situation 
is  more  complex,  as  Figs.  10  and  11  show. 

(c)  since  the  GISS  model  predicts  Ri{ cr)  «  0(1),  it  gives  rise  to  mixed  layer  depths  m  stable  sit¬ 
uations  in  agreement  with  the  data  (Burchard  and  Bolding,  2001;  Burchard  and  Deleersmjder, 
2001).  At  the  same  time,  in  unstable  convective-like  situations,  the  predictions  for  the  mixed 
layer  depths  are  closer  to  the  LDO  data  than  those  of  the  other  two  models.  In  so  doing, 
it  avoids  the  shortcomings  of  both  MY  and  KPP  models. 

The  GISS  model  is  still  missing  two  important  ingredients,  however.  The  first  is  rotation.  The 
Reynolds  stress  method  on  which  the  GISS  model  is  based,  has  definite  prescriptions  to  include 
the  effect  of  Q.  The  new  rotation-dependent  diffusivities  have  already  been  derived  (in  the  local 
limit)  and  they  will  be  used  next.  The  second  missing  ingredient  is  non-locality  (counter-gradients) 
represented  by  TOMs  like  the  velocity  skewness  that  characterizes  the  topology  of  the  descending 
and  ascending  plumes  and  for  which  measurements  are  available  to  test  the  model  results. 
Heuristic  models  cannot  compute  the  TOMs  (which  must  be  borrowed  from  other  sources)  while 
the  RSM  offers  a  well  defined  procedure.  Considerable  work  to  derive  analytic  expressions  for  the 
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TOMs  has  already  been  done  (Canuto  et  al.,  1994,  2001a, b,  submitted  for  publication)  and  their 
inclusion  in  a  DC  model  is  now  possible. 

Thirdly,  in  both  KPP  and  GISS  schemes,  lateral  advection  that  dominates  the  final  stages  of 
deep  convection  has  thus  far  been  modeled  using  the  heuristic  Gent-McWilliams  parameteriza¬ 
tion.  Recent  work  (Canuto  and  Dubovikov,  submitted  for  publication)  has  shown  that  a  meso- 
scale  model  can  be  derived  from  the  basic  dynamic  equations  (Navier-Stokes  and  scalar 
equations).  The  new  model  will  hopefully  improve  the  traditional  GM  model.  The  new  model  has 
however  not  yet  been  independently  tested.  When  that  is  done  and  if  significant  improvements 
over  the  GM  model  will  emerge,  the  lateral  advection  part  of  the  DC  model  will  also  be  improved. 
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Appendix  A.  Numerical  scheme 

The  OGCMs  were  integrated  for  126+  momentum  years.  In  Cl-2,  for  the  first  96  years,  a  depth 
acceleration  was  used  so  that  the  tracer  timestep  increased  from  10  times  the  momentum  timestep 
at  the  surface  to  100  times  the  momentum  timestep  in  the  abyss.  We  have  recently  found  that 
using  the  surface  tracer  timestep  throughout  the  column  for  the  first  96  years  makes  a  noticeable 
improvement  in  our  model,  in  which  the  deep  ocean  diffusivities  vary,  though  not  in  the  KPP 
model  with  constant  background  diffusivities.  In  Large  et  al.  (1997),  the  code  was  run  with  the  last 
30  years  of  the  simulations  using  a  tracer  timestep  equal  to  the  momentum  timestep.  However, 
Kamenkovich  et  al.  (2002)  have  shown  that  use  of  a  depth-independent  tracer  acceleration  factor 
of  12  did  not  degrade  results  vis-a-vis  use  of  equal  time-steps  in  the  MOM2  ocean  component  of 
their  climate  model;  we  have  found,  in  a  trial  with  the  NCOM,  negligible  difference  in  the  basin- 
averaged  temperature  and  salinity  profiles  between  a  run  in  which  the  last  30  years  used  a  depth- 
independent  tracer  acceleration  of  a  factor  of  10  and  one  in  which  the  last  30  years  were  run 
synchronously.  The  results  reported  here  for  both  the  KPP  and  GISS  models  accordingly  use  a 
tracer  timestep  10  times  the  momentum  timestep  throughout  the  ocean  for  the  entire  126-year  run. 
In  the  simulations  Cl-2,  an  artificial  cap  of  0.1  m2  s-1  was  placed  on  the  diffusivities  calculated  by 
the  GISS  model  due  to  fears  of  numerical  instability  that  turned  out  to  be  unfounded.  For  the 
runs  reported  here,  the  cap  was  raised  to  100  m2s_1,  a  limit  which  was  rarely  reached  in  the 
simulation  and  which  is  larger  than  the  diffusivities  predicted  for  the  Labrador  Sea.  With 
the  exceptions  of  the  changes  in  timestep  and  GISS  mixing  model,  the  OGCM  code  was  the  same 
as  in  C2. 
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